Code Packages:
# Load packages
library(gridExtra) # For using grid.table to save tables as images
## Warning: package 'gridExtra' was built under R version 4.4.3
library(flextable) # Alternative method to save df as image
## Warning: package 'flextable' was built under R version 4.4.3
library(tidyverse)
## Warning: package 'ggplot2' was built under R version 4.4.3
## Warning: package 'lubridate' was built under R version 4.4.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine() masks gridExtra::combine()
## ✖ purrr::compose() masks flextable::compose()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(zoo)
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(lme4)
## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
library(lmerTest)
##
## Attaching package: 'lmerTest'
##
## The following object is masked from 'package:lme4':
##
## lmer
##
## The following object is masked from 'package:stats':
##
## step
library(rlang)
##
## Attaching package: 'rlang'
##
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
library(nortest) # For ad.test
library(car) # For homogeneity of variance with leveneTest
## Warning: package 'car' was built under R version 4.4.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.4.3
##
## Attaching package: 'car'
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## some
library(rcompanion) # For scheirerRayHare test
## Warning: package 'rcompanion' was built under R version 4.4.3
library(pastecs) # For stat.desc
## Warning: package 'pastecs' was built under R version 4.4.3
##
## Attaching package: 'pastecs'
##
## The following objects are masked from 'package:dplyr':
##
## first, last
##
## The following object is masked from 'package:tidyr':
##
## extract
library(ggpubr) # For assumption plots
## Warning: package 'ggpubr' was built under R version 4.4.3
##
## Attaching package: 'ggpubr'
##
## The following objects are masked from 'package:flextable':
##
## border, font, rotate
library(conover.test) # For the conover.Iman() test function
library(DescTools) # Contain alternative ConoverTest function
## Warning: package 'DescTools' was built under R version 4.4.3
##
## Attaching package: 'DescTools'
##
## The following object is masked from 'package:car':
##
## Recode
library(rstatix) # tidyverse adapted stat tests and mshapiro_test() functions
## Warning: package 'rstatix' was built under R version 4.4.3
##
## Attaching package: 'rstatix'
##
## The following object is masked from 'package:stats':
##
## filter
library(DHARMa) # For multi-level lm models and shapiro test
## Warning: package 'DHARMa' was built under R version 4.4.3
## This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(patchwork) # For combining ggplots
## Warning: package 'patchwork' was built under R version 4.4.3
# Read and call data into df
# Primary samples df
ABL90 <- read.csv("Raw_Data/ABL90_Raw.csv")
# River's Parturition metadata
partuition_subcat <- read.csv("Raw_Data/River's_Rockfish_Metadata_Parturition_V7.csv")
# Data sifting: ABL90 dataset
# Step 1: Remove missing info
ABL_set1 <- ABL90 %>%
filter(Patient.ID_edited != "")
# Step 2: Separate Patient.ID into columns of sample types: blood plasma and instant freeze plasma
ABL_set2 <- separate(ABL_set1, 'Patient.ID_edited', into = c("Patient.ID_edited", "Sample_Type"), sep = ",")
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 4 rows [489, 619, 1034,
## 1177].
# Step 3:
# Convert Patient.ID_edited column data to character type
ABL_set2$Patient.ID_edited <- as.character(ABL_set2$Patient.ID_edited)
# Step 4 & 5: Sussing out specific sample errors
# Remove insufficient samples
ABL_set3 <- ABL_set2 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Insufficient sample"))
# Remove inhomogeneous samples
ABL_set3 <- ABL_set3 %>%
filter(!is.na("Type")) %>%
filter(!str_detect(Errors.detected.during.measurement, "Inhomogeneous sample"))
# Step 6: Filter for only blood samples
ABL_b_samp <- ABL_set3 %>%
filter(Sample_Type == "b")
Join ABL90 df with parturition_subcat df: Making ABL_merged
# Rename columns to match: Change 'ID' in parturition_subcat to 'Patient.ID_edited' so it is ready to join with ABL90 df
# Please note the 'Treatment' col in parturition_subcat does not match with the finalized 'Ambient.Or.OAH' col in metadata_atresia_guide or ABL90_merged df
# Rename 'ID' col to 'patient.ID_edited'
partuition_subcat <- partuition_subcat %>%
rename(Patient.ID_edited = ID)
# Connect main ABL90 dataset with my (River's) sub categorized parturition metadata
# ABL_merged <- partuition_subcat %>%
# left_join(ABL_b_samp, by = "Patient.ID_edited")
ABL_merged <- ABL_b_samp %>%
inner_join(partuition_subcat, by = "Patient.ID_edited")
Rename Columns Neatly:
# Rename Treatment & Parturition Outcome
ABL_merged <- ABL_merged %>%
rename(Ambient.Or.OAH = Consensus_General_Treatment,
Pregnant.Or.Atresia = Consensus_Atresia_Or_Pregnant)
Remove Replicate ID’s
# Cross Check all Columns: Cross validation to suss correct replicate row
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9782D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9782D 12/28/2023 10:04 1016 C 65uL Ambient
## 2 9782D 1/9/2024 17:21 863 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.348 0.5 9.8 184 149
## 2 Post Parturition 7.446 0.7 13.3 181 161
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 2.8 1.32 4.2
## 2 2.9 1.31 6.8
check_params <- ABL_merged %>%
select(Patient.ID_edited, Time, Sample.., Measuring.Mode, Ambient.Or.OAH, Pregnant.Or.Atresia, pH, Glu..mmol.L., Hct...., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L., pCO2..mmHg.) %>%
filter(Patient.ID_edited == "9783D")
print(check_params)
## Patient.ID_edited Time Sample.. Measuring.Mode Ambient.Or.OAH
## 1 9783D 1/19/2024 10:37 934 S 65uL Ambient
## 2 9783D 1/7/2024 12:32 860 S 65uL Ambient
## Pregnant.Or.Atresia pH Glu..mmol.L. Hct.... Na...mmol.L. Cl...mmol.L.
## 1 Post Parturition 7.484 1.7 0.2 175 193
## 2 Post Parturition 7.450 1.3 14.7 179 164
## K...mmol.L. Ca....mmol.L. pCO2..mmHg.
## 1 34.5 1.29 8.5
## 2 2.5 1.30 5.3
# Removed replicates: 9782D (2x) and 9783D (2x)
ABL_merged <- ABL_merged %>%
filter(Sample.. != "863",
Sample.. != "934")
Review Replicates:
| Row | ID | Time | Sample.. | Measuring.Mode | pH | Glu.mmol.L. | Hct…. | Na…mmol.L. | Cl…mmol.L. | K…mmol.L. | Ca.mmol.L. |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 32 | 9782D | 12/28/2023 10:04 | 1016 | C 65uL | 7.348 | 0.5 | 9.8 | 184 | 149 | 2.8 | 1.32 |
| 39 | 9782D | 1/9/2024 17:21 | 863 | S 65uL | 7.446 | 0.7 | 13.3 | 181 | 161 | 2.9 | 1.31 |
| 36 | 9783D | 1/19/2024 10:37 | 934 | S 65uL | 7.484 | 1.7 | 0.2 | 175 | 193 | 34.5 | 1.29 |
| 41 | 9783D | 1/7/2024 12:32 | 860 | S 65uL | 7.450 | 1.3 | 14.7 | 179 | 164 | 2.5 | 1.30 |
Methods for replicate sample removal: Check processing date of IDs, and remove samples that do not match that date.
Case ID_9782D:
Origianl ABL90 set contains:
9782D,b (Sample.. 1016, Time = 12/28/2023 10:04:00 AM)
9782D,p (Sample.. 1020, Time = 12/28/2023 10:25:00 AM)
9782D,p (Sample.. 865, Time = 1/9/2024 5:29:00 PM)
9782D,b (Sample.. 863, Time = 1/9/2024 5:21:00 PM)
2024_Fish_Metadata.csv says ID_9782D was processed 3/29/2024
Note Processing date in metadata will not match time date stamp in ABL90 data due to the machine being reset, and innacurrate dates being automatically recorded.
Deciding info: - For ID-9782D, refer to metadata sheet and look for other samples with that same processing date, then view thoes samples ABL90 ‘Time’. Try to find a matching ‘Time’ with other IDs and 9782D that were processed together during the same time recorded by the machine.
ID_97838,b (single obs, Sample.. 1024) Time = 12/28/2023 11:49:00 AM
Since 9782D,b Sample.. 1016 and 97838,b Sample.. 1024 both have a date timestamp occurrence on 12/28/2023, I have decided to keep 9782D,b Sample.. 1016 and drop 9783D,b Sample.. 863
Case ID_9783D
Original ABL90 dataset contains:
9783D,b (Sample.. 934, Time = 1/19/2024 10:37:00 AM)
9783D,c (Sample.. 861, Time = 1/7/2024 12:38:00 PM)
9783D,b (Sample.. 860, Time = 1/7/2024 12:32:00 PM)
2024_Fish_Metadata.csv says ID_9783D was processed 3/20/2024
Informed that ID_9783D: Sample.. 934 is a plasma sample (on account of its low hct value), remove this and keep Sample.. 860 (parameter values appear correct).
Remove Motalities and No info IDs and assign analysis ready data-frames:
# New df with Moralities removed: Note none of these samples made it into ABL90 df anyway, so looks like they are already filtered out.
Live_Samples <- ABL_merged %>%
filter(Patient.ID_edited != "9780C", # Mortality
Patient.ID_edited != "777AE", # Mortality
Patient.ID_edited != "777CA") # Mortality after parturition
# New df with mortality and 'No info' Id's removed
Primary_Samples <- Live_Samples %>%
filter(Patient.ID_edited != "777A0", # No info
Patient.ID_edited != "9782F", # No info
Patient.ID_edited != "777B3", # No info
Patient.ID_edited != "777AA", # No info
Patient.ID_edited != "777DE", # No info
Patient.ID_edited != "777CE", # No info
Patient.ID_edited != "777A6") # No info
# New df of Only Ambient Treatment: For testing parturition success
Ambient_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient")
# New df of Fecundity Samples:
Fecundity_Samples <- Primary_Samples %>%
filter(Fecundity_Or_Brood_Count != "NA",
Fecundity_Class != "NA") # removes 97706
# Fecundity samples without atresia Ids
Fecundity_No_Atresia_Samples <- Primary_Samples %>%
filter(All_Fecundity != "Na",
All_Fecundity != 0) # removes atresia IDs
# Fecundity Ambient Samples
Amb_Fec_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient",
All_Fecundity != "NA") # removes OAH treatment samples and ID 97706
# Fecundity Ambient Samples without atresia Ids
Amb_Fec_No_Atresia_Samples <- Primary_Samples %>%
filter(Ambient.Or.OAH == "Ambient", # removes OAH treatment
All_Fecundity != "NA", # removes all missing info IDs
All_Fecundity != 0) # removes atresia IDs
Change Data Types:
# Change Continuous Data to type to numeric, double, or factor
# Change data type of ions to factor or numeric
#glimpse(General_Samples)
#glimpse(Ambient_Samples)
Primary_Samples <- Primary_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Ambient_Samples <- Ambient_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Fecundity_No_Atresia_Samples <- Fecundity_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_Samples <- Amb_Fec_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
mutate(across(c(All_Fecundity, pH, Hct...., Glu..mmol.L., Na...mmol.L., Cl...mmol.L., K...mmol.L., Ca....mmol.L.), as.numeric))
# glimpse(Primary_Samples)
# glimpse(Ambient_Samples)
# glimpse(Fecundity_Samples)
# glimpse(Amb_Fec_Samples)
# Fecundity data
#unique(Fecundity_Samples$Fecundity_Class)
# Arrange the order of parturition categories for visualizations & tests
# Primary Samples
# Arrange Treatment
Primary_Samples$Ambient.Or.OAH <- ordered(Primary_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
# Arrange Parturition outcome
Primary_Samples$Pregnant.Or.Atresia <- ordered(Primary_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
# Arrange Brood Condition
Primary_Samples$Consensus_Brood_Condition <- ordered(Primary_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Very Low", "Atresia"))
# Arrange Fecundity Class
Primary_Samples$Fecundity_Class <- ordered(Primary_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Primary_Samples <- Primary_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Ambient data
Ambient_Samples$Ambient.Or.OAH <- ordered(Ambient_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
Ambient_Samples$Pregnant.Or.Atresia <- ordered(Ambient_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
Ambient_Samples$Consensus_Brood_Condition <- ordered(Ambient_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))
Ambient_Samples$Fecundity_Class <- ordered(Ambient_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Ambient_Samples <- Ambient_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Ambient Fecundity data
Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH <- ordered(Amb_Fec_No_Atresia_Samples$Ambient.Or.OAH, levels = c("Ambient", "OAH"))
Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia <- ordered(Amb_Fec_No_Atresia_Samples$Pregnant.Or.Atresia, levels = c("Post Parturition", "Atresia"))
Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition <- ordered(Amb_Fec_No_Atresia_Samples$Consensus_Brood_Condition, levels = c("Excellent", "Normal", "Low", "Atresia"))
Amb_Fec_No_Atresia_Samples$Fecundity_Class <- ordered(Amb_Fec_No_Atresia_Samples$Fecundity_Class, levels = c("High (>50,000)", "Low (~1,000s)", "Atresia"))
# Remove NAs from ordered character factors
Amb_Fec_No_Atresia_Samples <- Amb_Fec_No_Atresia_Samples %>%
filter(Fecundity_Class != is.na(Fecundity_Class)) %>%
filter(Consensus_Brood_Condition != is.na(Consensus_Brood_Condition))
# Save merged dataset in a worked folder
write.csv(x = Ambient_Samples, file = "Worked-Data/Ambient_Samples")
# Sample Size
# Ambient data:
ambient_sample_size_table <- Ambient_Samples %>%
group_by(Ambient.Or.OAH, Consensus_Brood_Condition) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_sample_size_table)
## # A tibble: 4 × 3
## Ambient.Or.OAH Consensus_Brood_Condition n_size
## <ord> <ord> <int>
## 1 Ambient Excellent 1
## 2 Ambient Normal 4
## 3 Ambient Low 3
## 4 Ambient Atresia 13
ambient_fecundity_class_sample_size_table <- Ambient_Samples %>%
group_by(Ambient.Or.OAH, Fecundity_Class) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(ambient_fecundity_class_sample_size_table)
## # A tibble: 3 × 3
## Ambient.Or.OAH Fecundity_Class n_size
## <ord> <ord> <int>
## 1 Ambient High (>50,000) 4
## 2 Ambient Low (~1,000s) 4
## 3 Ambient Atresia 13
amb_fec_no_atresia_sample_size_table <- Amb_Fec_No_Atresia_Samples %>%
group_by(Ambient.Or.OAH, Weight_Adjusted_Fecundity_Fecundity.Mean_Weight) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop")
print(amb_fec_no_atresia_sample_size_table)
## # A tibble: 8 × 3
## Ambient.Or.OAH Weight_Adjusted_Fecundity_Fecundity.Mean_Weight n_size
## <ord> <dbl> <int>
## 1 Ambient 6.22 1
## 2 Ambient 68.8 1
## 3 Ambient 114. 1
## 4 Ambient 114. 1
## 5 Ambient 148. 1
## 6 Ambient 161. 1
## 7 Ambient 204. 1
## 8 Ambient 356. 1
Sample Size Barplot:
# View Sample Distributions
# Sample Size barplot
# Ambient Samples
# Consensus brood condition sample barplot
ambient_sample_size_barplot <- Ambient_Samples %>%
group_by(Consensus_Brood_Condition, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Consensus_Brood_Condition, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
labs(title = "Ambient Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot)
# Fecundity class sample barplot
ambient_sample_size_barplot2 <- Ambient_Samples %>%
group_by(Fecundity_Class, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Fecundity_Class, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(0, 3, 6, 9, 12, 15)) +
labs(title = "Ambient Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(ambient_sample_size_barplot2)
# Amb Fec No atresia sample size
amb_fec_no_atresia_barplot <- Amb_Fec_No_Atresia_Samples %>%
group_by(Fecundity_Class, Ambient.Or.OAH) %>%
summarise(n_size = n_distinct(Patient.ID_edited), .groups = "drop") %>%
ggplot(aes(x = Fecundity_Class, y = n_size)) +
geom_bar(stat = "identity") +
facet_grid(~ Ambient.Or.OAH) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 15.5), breaks = c(4, 8, 12)) +
labs(title = "Ambient No Atresia Sample Size",
x = "Parturition Outcome",
y = "Sample Size (n)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(amb_fec_no_atresia_barplot)
CO2 Boxplot
pCO2_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pCO2..mmHg.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = pCO2..mmHg.), alpha = 0.5, colour = "black") +
labs(title = "Carbon Dioxide", x = "Fecundity Class", y = "pCO2 (mmHg)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(pCO2_ambient_fec_class_boxplot)
# Ambient No Atresia Samples Scatterplot
pCO2_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pCO2..mmHg.), shape = 21, fill = "#189bff") +
labs(title = "Carbon Dioxide", x = "Weight Adjusted Fecundity", y = "pCO2 (mmHg)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pCO2_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# pCO2
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
pCO2_aov_ambient_fec <- aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
summary(pCO2_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 10.10 5.050 1.817 0.191
## Residuals 18 50.02 2.779
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(pCO2_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -1.2750000 -4.2833026 1.733303 0.5371486
## Atresia-High (>50,000) 0.5384615 -1.8940746 2.970998 0.8401814
## Atresia-Low (~1,000s) 1.8134615 -0.6190746 4.245998 0.1666659
# Install flexplot
# devtools::install_github("dustinfife/flexplot", ref="development")
require(flexplot)
## Loading required package: flexplot
##
## Attaching package: 'flexplot'
## The following object is masked from 'package:ggplot2':
##
## flip_data
pCO2_ambient_fec_sim <- simulateResiduals(fittedModel = pCO2_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(pCO2_ambient_fec_sim)
plotResiduals(pCO2_ambient_fec_sim)
testDispersion(pCO2_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.92221, p-value = 0.888
## alternative hypothesis: two.sided
testOutliers(pCO2_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: pCO2_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(pCO2_aov_ambient_fec)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 7 2.080913 0.052875 NA
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$pCO2..mmHg.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$pCO2..mmHg.
## W = 0.97964, p-value = 0.9201
shapiro.test(residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.95328, p-value = 0.3922
# QQplot:
pCO2_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pCO2 Residual QQplot",
subtitle = "residuals(aov(pCO2 ~ Fecundity_Class, data = Ambient_Samples))",
x = "pCO2 Theoretical Quantiles (Predicted values)",
y = "pCO2 Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(pCO2_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: pCO2..mmHg. by Fecundity_Class
## Bartlett's K-squared = 0.88695, df = 2, p-value = 0.6418
# Nonparametric variance test (Levene's test):
leveneTest(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5105 0.6086
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: pCO2..mmHg. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.6912, df = 2, p-value = 0.2604
# Nonparametric Post Hoc: Dunn's test
DunnTest(pCO2..mmHg. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) -4.125000 0.6921
## Atresia-High (>50,000) 1.673077 0.6921
## Atresia-Low (~1,000s) 5.798077 0.3043
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# pCO2
# Linear Model
pCO2_amb_fec_no_atresia_lm <- lm(pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(pCO2_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.760 -1.069 0.120 1.401 1.880
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 4.770067 1.138478 4.190
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.001314 0.006485 0.203
## Pr(>|t|)
## (Intercept) 0.00575 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.84608
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.777 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.006801, Adjusted R-squared: -0.1587
## F-statistic: 0.04108 on 1 and 6 DF, p-value: 0.8461
pCO2_ambient_sim <- simulateResiduals(fittedModel = pCO2_amb_fec_no_atresia_lm)
plot(pCO2_ambient_sim)
plotResiduals(pCO2_ambient_sim)
testDispersion(pCO2_ambient_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(pCO2_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: pCO2_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(pCO2_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 21 -2.262971 0.073079 0.58463
# Normality test
shapiro.test(residuals(pCO2_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(pCO2_amb_fec_no_atresia_lm)
## W = 0.93831, p-value = 0.5945
# qq plot
ggqqplot(residuals(lm(pCO2..mmHg. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
Residuals Plot (if significant)
# pCO2
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# Original lm model (no transformation)
pCO2_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pCO2..mmHg.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Carbon Dioxide", x = "Weight Adjusted Fecundity", y = "pCO2 (mmHg)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pCO2_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
pH boxplot
pH_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pH), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Fecundity Class", y = "pH") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(pH_ambient_fec_class_boxplot)
# Ambient No Atresia Samples Scatterplot
pH_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH), shape = 21, fill = "#189bff") +
labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# pH
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
pH_aov_ambient_fec <- aov(pH ~ Fecundity_Class, data = Ambient_Samples)
summary(pH_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.02183 0.010916 1.336 0.288
## Residuals 18 0.14711 0.008173
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(pH_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = pH ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.02675000 -0.1363959 0.18989587 0.9084726
## Atresia-High (>50,000) -0.05080769 -0.1827287 0.08111329 0.5966404
## Atresia-Low (~1,000s) -0.07755769 -0.2094787 0.05436329 0.3141611
pH_ambient_fec_sim <- simulateResiduals(fittedModel = pH_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(pH_ambient_fec_sim)
plotResiduals(pH_ambient_fec_sim)
testDispersion(pH_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(pH_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: pH_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(pH_aov_ambient_fec)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 13 3.11884 0.0062471 0.13119
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$pH)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$pH
## W = 0.98534, p-value = 0.9808
shapiro.test(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97088, p-value = 0.7524
# QQplot:
pH_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient pH Residual QQplot",
subtitle = "residuals(aov(pH ~ Fecundity_Class, data = Ambient_Samples))",
x = "pH Theoretical Quantiles (Predicted values)",
y = "pH Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(pH_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: pH by Fecundity_Class
## Bartlett's K-squared = 1.4521, df = 2, p-value = 0.4838
# Nonparametric variance test (Levene's test):
leveneTest(pH ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.44 0.6508
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: pH by Fecundity_Class
## Kruskal-Wallis chi-squared = 3.0375, df = 2, p-value = 0.219
# Nonparametric Post Hoc: Dunn's test
DunnTest(pH ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) 2.250000 0.6424
## Atresia-High (>50,000) -3.519231 0.6424
## Atresia-Low (~1,000s) -5.769231 0.3117
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Ph
# Linear Model
pH_amb_fec_no_atresia_lm <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(pH_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.11361 -0.02694 0.01591 0.04102 0.05887
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 7.4734874 0.0403735 185.109
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -0.0002228 0.0002300 -0.969
## Pr(>|t|)
## (Intercept) 1.68e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.37
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06302 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.1352, Adjusted R-squared: -0.008899
## F-statistic: 0.9383 on 1 and 6 DF, p-value: 0.3701
pH_ambient_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm)
plot(pH_ambient_sim)
## qu = 0.25, log(sigma) = -2.642568 : outer Newton did not converge fully.
plotResiduals(pH_amb_fec_no_atresia_lm)
## qu = 0.25, log(sigma) = -2.642568 : outer Newton did not converge fully.
testDispersion(pH_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(pH_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: pH_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(pH_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 3 -2.862653 0.035297 0.28237
# Normality test
shapiro.test(residuals(pH_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(pH_amb_fec_no_atresia_lm)
## W = 0.90496, p-value = 0.32
# qq plot
ggqqplot(residuals(lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root (undo log scale?)
pH_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(pH = sqrt(pH))
# Rerun lm model with square root transformed data
pH_amb_fec_no_atresia_lm_sqrt <- lm(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = pH_ambient_sqrt)
summary(pH_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = pH_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.020889 -0.004920 0.002934 0.007526 0.010787
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.734e+00 7.412e-03 368.821
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -4.078e-05 4.222e-05 -0.966
## Pr(>|t|)
## (Intercept) 2.68e-14 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.371
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01157 on 6 degrees of freedom
## Multiple R-squared: 0.1346, Adjusted R-squared: -0.009642
## F-statistic: 0.9331 on 1 and 6 DF, p-value: 0.3714
pH_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = pH_amb_fec_no_atresia_lm_sqrt)
plot(pH_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.
plotResiduals(pH_amb_fec_no_atresia_lm_sqrt)
## qu = 0.25, log(sigma) = -2.642951 : outer Newton did not converge fully.
testDispersion(pH_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
pH_amb_fec_no_atresia_lm_res <- residuals(pH_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(pH_amb_fec_no_atresia_lm_res) # normal
##
## Shapiro-Wilk normality test
##
## data: pH_amb_fec_no_atresia_lm_res
## W = 0.90414, p-value = 0.3147
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(pH ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Fecundity_No_Atresia_Samples)
Residuals Plot (if significant)
# pH
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
# Original lm model (no transformation)
pH_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Weight Adjusted Fecundity", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(pH_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
# hct
# Hct boxplot
hct_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Hct....), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Fecundity", y = "Hematocrit (%)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(hct_ambient_fec_class_boxplot)
# Ambient No Atresia Sample Scatterplot
hct_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....), shape = 21, fill = "#189bff") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# hct
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
hct_aov_ambient_fec <- aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
summary(hct_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 73.83 36.92 4.824 0.021 *
## Residuals 18 137.74 7.65
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(hct_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -3.4 -8.3921454 1.592145 0.2186627
## Atresia-High (>50,000) 1.5 -2.5366864 5.536686 0.6176987
## Atresia-Low (~1,000s) 4.9 0.8633136 8.936686 0.0162654
hct_ambient_fec_sim <- simulateResiduals(fittedModel = hct_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(hct_ambient_fec_sim)
plotResiduals(hct_ambient_fec_sim)
testDispersion(hct_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(hct_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: hct_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(hct_aov_ambient_fec)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 20 2.329698 0.032404 0.68048
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Hct....)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Hct....
## W = 0.97482, p-value = 0.8355
shapiro.test(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.97094, p-value = 0.7538
# QQplot:
hct_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Hematocrit Residual QQplot",
subtitle = "residuals(aov(Hct.... ~ Fecundity_Class, data = Ambient_Samples))",
x = "Hct Theoretical Quantiles (Predicted values)",
y = "Hct Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(hct_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Hct.... by Fecundity_Class
## Bartlett's K-squared = 1.7219, df = 2, p-value = 0.4228
# Nonparametric variance test (Levene's test):
leveneTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2003 0.3241
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Hct.... by Fecundity_Class
## Kruskal-Wallis chi-squared = 6.7821, df = 2, p-value = 0.03367
# Nonparametric Post Hoc: Dunn's Test
DunnTest(Hct.... ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) -6.750000 0.2476
## Atresia-High (>50,000) 2.480769 0.4843
## Atresia-Low (~1,000s) 9.230769 0.0277 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# hct
# Linear Model
hct_amb_fec_no_atresia_lm <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(hct_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.5587 -1.5142 -0.7795 0.5244 7.1866
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 12.81211 2.32735 5.505
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.01494 0.01326 1.127
## Pr(>|t|)
## (Intercept) 0.00151 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.30267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.633 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.1748, Adjusted R-squared: 0.03725
## F-statistic: 1.271 on 1 and 6 DF, p-value: 0.3027
hct_ambient_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm)
plot(hct_ambient_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## We had to increase `err` for some of the quantiles. See fit$calibr$err
plotResiduals(hct_amb_fec_no_atresia_lm)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## We had to increase `err` for some of the quantiles. See fit$calibr$err
testDispersion(hct_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(hct_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: hct_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(hct_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 20 3.953321 0.010815 0.086518
# Normality test
shapiro.test(residuals(hct_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(hct_amb_fec_no_atresia_lm)
## W = 0.86076, p-value = 0.1222
# qq plot
ggqqplot(residuals(lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root (undo log scale?)
hct_ambient_sqrt <- Ambient_Samples %>%
mutate(Hct.... = sqrt(Hct....))
# Rerun lm model with square root transformed data
hct_amb_fec_no_atresia_lm_sqrt <- lm(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
summary(hct_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = hct_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45751 -0.18537 -0.07738 0.06130 0.87682
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 3.549914 0.292111 12.153
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.002034 0.001664 1.223
## Pr(>|t|)
## (Intercept) 1.89e-05 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.456 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.1995, Adjusted R-squared: 0.06603
## F-statistic: 1.495 on 1 and 6 DF, p-value: 0.2673
hct_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = hct_amb_fec_no_atresia_lm_sqrt)
plot(hct_amb_fec_sqrt_sim)
plotResiduals(hct_amb_fec_no_atresia_lm_sqrt)
testDispersion(hct_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
hct_amb_sqrt_res <- residuals(hct_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(hct_amb_sqrt_res)
##
## Shapiro-Wilk normality test
##
## data: hct_amb_sqrt_res
## W = 0.87306, p-value = 0.1615
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Hct.... ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = hct_ambient_sqrt)
Residuals Plot (if significant)
# hct
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
hct_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecundity", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(hct_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Glucose graphs
# Glucose
# Glucose boxplot
glucose_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(glucose_ambient_fec_class_boxplot)
# Glucose Ambient No Atresia Scatterplot
glucose_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.), shape = 21, fill = "#189bff") +
labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# glucose
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
glucose_aov_ambient_fec <- aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(glucose_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 1.208 0.6041 1.317 0.293
## Residuals 18 8.257 0.4587
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(glucose_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.125 -1.0973103 1.347310 0.9632210
## Atresia-High (>50,000) 0.550 -0.4383693 1.538369 0.3518867
## Atresia-Low (~1,000s) 0.425 -0.5633693 1.413369 0.5277964
glu_ambient_fec_sim <- simulateResiduals(fittedModel = glucose_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(glu_ambient_fec_sim)
plotResiduals(glu_ambient_fec_sim)
testDispersion(glu_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(glu_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: glu_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(glucose_aov_ambient_fec)
## rstudent unadjusted p-value Bonferroni p
## 13 4.830979 0.0001563 0.0032822
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Glu..mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Glu..mmol.L.
## W = 0.85634, p-value = 0.005473
shapiro.test(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.82056, p-value = 0.001376
# QQplot:
glucose_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Glucose Interactive Residual QQplot",
subtitle = "residuals(aov(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Glucose Theoretical Quantiles (Predicted values)",
y = "Glucose Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(glucose_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Glu..mmol.L. by Fecundity_Class
## Bartlett's K-squared = 3.669, df = 2, p-value = 0.1597
# Nonparametric variance test (Levene's test):
leveneTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.9391 0.4093
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Glu..mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 4.0755, df = 2, p-value = 0.1303
# Nonparametric Post Hoc: Dunn's test
DunnTest(Glu..mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) 2.625000 0.5480
## Atresia-High (>50,000) 6.663462 0.1778
## Atresia-Low (~1,000s) 4.038462 0.5060
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# glucose
# Linear Model
glucose_amb_fec_no_atresia_lm <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(glucose_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71290 -0.36526 -0.05942 0.16524 1.19856
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.5622030 0.4085226 3.824
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0003435 0.0023269 0.148
## Pr(>|t|)
## (Intercept) 0.00872 **
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.88746
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6377 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.00362, Adjusted R-squared: -0.1624
## F-statistic: 0.0218 on 1 and 6 DF, p-value: 0.8875
glucose_ambient_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm)
plot(glucose_ambient_sim)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
plotResiduals(glucose_amb_fec_no_atresia_lm)
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
testDispersion(glucose_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(glucose_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: glucose_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(glucose_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 5 3.287331 0.021778 0.17422
# Normality test:
shapiro.test(residuals(glucose_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(glucose_amb_fec_no_atresia_lm)
## W = 0.92103, p-value = 0.4383
# qq plot
ggqqplot(residuals(lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root
glucose_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Glu..mmol.L. = sqrt(Glu..mmol.L.))
# Rerun lm model with square root transformed data
glucose_amb_fec_no_atresia_lm_sqrt <- lm(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
summary(glucose_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = glucose_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30380 -0.13647 -0.01130 0.07862 0.42723
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.2241905 0.1551624 7.890
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001917 0.0008838 0.217
## Pr(>|t|)
## (Intercept) 0.00022 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.83544
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2422 on 6 degrees of freedom
## Multiple R-squared: 0.007783, Adjusted R-squared: -0.1576
## F-statistic: 0.04706 on 1 and 6 DF, p-value: 0.8354
glucose_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = glucose_amb_fec_no_atresia_lm_sqrt)
plot(glucose_amb_fec_sqrt_sim)
plotResiduals(glucose_amb_fec_no_atresia_lm_sqrt)
testDispersion(glucose_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
glucose_amb_fec_no_atresia_lm_res <- residuals(glucose_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(glucose_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: glucose_amb_fec_no_atresia_lm_res
## W = 0.96004, p-value = 0.8104
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Glu..mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = glucose_ambient_sqrt)
Residuals Plot (if significant)
# glucose
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
glucose_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Weight Adjusted Fecundity", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(glucose_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
# Na+
# sodium boxplot
sodium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(sodium_ambient_fec_class_boxplot)
# sodium ambient no atresia scatterplot
sodium_amb_fec_no_atresia_scatterplot <- Ambient_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), shape = 21, fill = "#189bff") +
labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_amb_fec_no_atresia_scatterplot)
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Differences: Between Fecundity Class
# Na+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
sodium_aov_ambient_fec <- aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(sodium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 156.8 78.38 0.865 0.438
## Residuals 18 1631.8 90.66
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(sodium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -2.000000 -19.18271 15.18271 0.9526465
## Atresia-High (>50,000) 4.480769 -9.41330 18.37484 0.6939179
## Atresia-Low (~1,000s) 6.480769 -7.41330 20.37484 0.4737484
na_ambient_fec_sim <- simulateResiduals(fittedModel = sodium_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(na_ambient_fec_sim)
plotResiduals(na_ambient_fec_sim)
testDispersion(na_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(na_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: na_ambient_fec_sim
## outliers at both margin(s) = 1, observations = 21, p-value = 0.1546
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.001204883 0.238159910
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0.04761905
outlierTest(sodium_aov_ambient_fec)
## rstudent unadjusted p-value Bonferroni p
## 17 12.20437 7.7681e-10 1.6313e-08
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Na...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Na...mmol.L.
## W = 0.52676, p-value = 3.625e-07
shapiro.test(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.56321, p-value = 8.235e-07
# QQplot:
sodium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Sodium Residual QQplot",
subtitle = "residuals(aov(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Sodium Theoretical Quantiles (Predicted values)",
y = "Sodium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(sodium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Na...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 7.4634, df = 2, p-value = 0.02395
# Nonparametric variance test (Levene's test):
leveneTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.2098 0.8127
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Na...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 2.8608, df = 2, p-value = 0.2392
# Nonparametric Post Hoc: Dunn's test
DunnTest(Na...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) -3.375000 0.8730
## Atresia-High (>50,000) 2.451923 0.8730
## Atresia-Low (~1,000s) 5.826923 0.2899
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Na+
# Linear Model
sodium_amb_fec_no_atresia_lm <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(sodium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7526 -1.0177 0.5473 2.5230 3.4276
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.754e+02 2.355e+00 74.474
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 2.289e-03 1.342e-02 0.171
## Pr(>|t|)
## (Intercept) 3.94e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.87
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.677 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.004827, Adjusted R-squared: -0.161
## F-statistic: 0.0291 on 1 and 6 DF, p-value: 0.8701
sodium_ambient_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm)
plot(sodium_ambient_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
plotResiduals(sodium_ambient_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
testDispersion(sodium_ambient_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(sodium_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: sodium_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(sodium_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 2 -2.997851 0.030174 0.24139
# Normailty test
shapiro.test(residuals(sodium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(sodium_amb_fec_no_atresia_lm)
## W = 0.89794, p-value = 0.2768
# qq plot
ggqqplot(residuals(lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root
sodium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Na...mmol.L. = sqrt(Na...mmol.L.))
# Rerun lm model with square root transformed data
sodium_amb_fec_no_atresia_lm_sqrt <- lm(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
summary(sodium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = sodium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25662 -0.03796 0.02115 0.09537 0.12926
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.324e+01 8.923e-02 148.42
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 8.630e-05 5.083e-04 0.17
## Pr(>|t|)
## (Intercept) 6.31e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.871
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1393 on 6 degrees of freedom
## Multiple R-squared: 0.004782, Adjusted R-squared: -0.1611
## F-statistic: 0.02883 on 1 and 6 DF, p-value: 0.8708
sodium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = sodium_amb_fec_no_atresia_lm_sqrt)
plot(sodium_amb_fec_sqrt_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
plotResiduals(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.75. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
testDispersion(sodium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
sodium_amb_fec_no_atresia_lm_res <- residuals(sodium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(sodium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: sodium_amb_fec_no_atresia_lm_res
## W = 0.89593, p-value = 0.2654
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Na...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = sodium_ambient_sqrt)
Residuals Plot (if significant)
# Na+
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
sodium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Weight Adjusted Fecundity", y = "Sodium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(sodium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Chloride graphs
# Cl-
# Chloride boxplot
chloride_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(chloride_ambient_fec_class_boxplot)
# Chloride Ambient No Atresia Scatterplot
chloride_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.), shape = 21, fill = "#189bff") +
labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# Cl-
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
chloride_aov_ambient_fec <- aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(chloride_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 3.0 1.504 0.084 0.92
## Residuals 18 322.2 17.902
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(chloride_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -0.5000000 -8.135556 7.135556 0.9847331
## Atresia-High (>50,000) 0.4615385 -5.712630 6.635707 0.9801564
## Atresia-Low (~1,000s) 0.9615385 -5.212630 7.135707 0.9170009
cl_ambient_fec_sim <- simulateResiduals(fittedModel = chloride_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(cl_ambient_fec_sim)
plotResiduals(cl_ambient_fec_sim)
testDispersion(cl_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(cl_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: cl_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(chloride_aov_ambient_fec)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 16 2.003759 0.061299 NA
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Cl...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Cl...mmol.L.
## W = 0.96502, p-value = 0.6224
shapiro.test(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.96555, p-value = 0.6339
# QQplot:
chloride_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Chloride Residual QQplot",
subtitle = "residuals(aov(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Chloride Theoretical Quantiles (Predicted values)",
y = "Chloride Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(chloride_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Cl...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 0.82247, df = 2, p-value = 0.6628
# Nonparametric variance test (Levene's test):
leveneTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.5546 0.5838
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Cl...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 0.4005, df = 2, p-value = 0.8185
# Nonparametric Post Hoc: Dunns
DunnTest(Cl...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) -2.0000000 1.0000
## Atresia-High (>50,000) 0.2115385 1.0000
## Atresia-Low (~1,000s) 2.2115385 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Remove Outlier
# Cl-
# Test outlier
cl_ambient_fec_sim <- simulateResiduals(fittedModel = chloride_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(cl_ambient_fec_sim)
plotResiduals(cl_ambient_fec_sim)
testDispersion(cl_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(cl_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: cl_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(chloride_aov_ambient_fec)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 16 2.003759 0.061299 NA
# Find Outlier
boxplot(Ambient_Samples$Cl...mmol.L., horizontal = TRUE)
ggplot(data = Ambient_Samples, aes(x = Fecundity_Class, y = Cl...mmol.L., fill = Ambient.Or.OAH)) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
geom_point(aes(color = Ambient.Or.OAH), position = position_dodge(width = 0.75), alpha = 0.3, color = "black")
# Remove Outlier
# cl_amb_no_outlier <- Ambient_Samples %>%
# filter(Cl...mmol.L. !=)
# Rerun Stat Tests
# chloride_aov_outlier_removed <- aov(Cl...mmol.L. ~ Fecundity_Class, data = cl_amb_no_outlier)
# summary(chloride_aov_outlier_removed)
#
#
# TukeyHSD(chloride_aov_outlier_removed)
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Cl-
# Linear Model
chloride_amb_fec_no_atresia_lm <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(chloride_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.7216 -2.2032 -0.5783 1.8178 5.0996
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 158.033787 2.144241 73.702
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -0.001938 0.012213 -0.159
## Pr(>|t|)
## (Intercept) 4.2e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.879
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.347 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.00418, Adjusted R-squared: -0.1618
## F-statistic: 0.02519 on 1 and 6 DF, p-value: 0.8791
chloride_ambient_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm)
plot(chloride_ambient_sim)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
plotResiduals(chloride_ambient_sim)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
testDispersion(chloride_ambient_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(chloride_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: chloride_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(chloride_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 21 2.177688 0.08135 0.6508
# Shapiro Test
shapiro.test(residuals(chloride_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(chloride_amb_fec_no_atresia_lm)
## W = 0.93939, p-value = 0.6051
# Residual QQ plot
ggqqplot(residuals(lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root
chloride_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Cl...mmol.L. = sqrt(Cl...mmol.L.))
# Rerun lm model with square root transformed data
chloride_amb_fec_no_atresia_lm_sqrt <- lm(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
summary(chloride_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = chloride_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.14856 -0.08721 -0.02279 0.07258 0.20200
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.257e+01 8.518e-02 147.566
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight -7.493e-05 4.852e-04 -0.154
## Pr(>|t|)
## (Intercept) 6.53e-12 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.882
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.133 on 6 degrees of freedom
## Multiple R-squared: 0.003959, Adjusted R-squared: -0.162
## F-statistic: 0.02385 on 1 and 6 DF, p-value: 0.8823
chloride_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = chloride_amb_fec_no_atresia_lm_sqrt)
plot(chloride_amb_fec_sqrt_sim)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
plotResiduals(chloride_amb_fec_no_atresia_lm_sqrt)
## qu = 0.75, log(sigma) = -2.939317 : outer Newton did not converge fully.
testDispersion(chloride_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(chloride_amb_fec_sqrt_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: chloride_amb_fec_sqrt_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(chloride_amb_fec_no_atresia_lm_sqrt)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 8 2.165434 0.082618 0.66094
# Residuals of model for shapiro test
chloride_amb_fec_no_atresia_lm_res <- residuals(chloride_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(chloride_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: chloride_amb_fec_no_atresia_lm_res
## W = 0.94043, p-value = 0.6153
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
#bartlett.test(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Cl...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = chloride_ambient_sqrt)
Residuals Plot (if significant)
# Cl-
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
chloride_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Weight Adjusted Fecundity", y = "Chloride (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(chloride_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
K+ graphs
# K+ Boxplot
potassium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(potassium_ambient_fec_class_boxplot)
# K+ Ambient No Atresia Scatterplot
potassium_amb_fec_no_atresia_scatterplot <- Amb_Fec_No_Atresia_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.), shape = 21, fill = "#189bff") +
labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_amb_fec_no_atresia_scatterplot)
Differences: Between Fecundity Class
# K+
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
potassium_aov_ambient_fec <- aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(potassium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.1113 0.05563 0.796 0.466
## Residuals 18 1.2583 0.06990
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(potassium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) 0.0250000 -0.4521380 0.5021380 0.9901955
## Atresia-High (>50,000) 0.1615385 -0.2242789 0.5473558 0.5449631
## Atresia-Low (~1,000s) 0.1365385 -0.2492789 0.5223558 0.6452587
k_ambient_fec_sim <- simulateResiduals(fittedModel = potassium_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(k_ambient_fec_sim)
plotResiduals(k_ambient_fec_sim)
testDispersion(k_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(k_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: k_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(potassium_aov_ambient_fec)
## rstudent unadjusted p-value Bonferroni p
## 13 3.878893 0.0012057 0.025319
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$K...mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$K...mmol.L.
## W = 0.81347, p-value = 0.001062
shapiro.test(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.88368, p-value = 0.0171
# QQplot:
potassium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Potassium Residual QQplot",
subtitle = "residuals(aov(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Potassium Theoretical Quantiles (Predicted values)",
y = "Potassium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(potassium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: K...mmol.L. by Fecundity_Class
## Bartlett's K-squared = 4.6368, df = 2, p-value = 0.09843
# Nonparametric variance test (Levene's test):
leveneTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 0.7708 0.4773
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric stat test: Kruskall Wallis
kruskal.test(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: K...mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 1.2081, df = 2, p-value = 0.5466
# Nonparametric Post Hoc: Dunn's test
DunnTest(K...mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) 0.375000 1.0000
## Atresia-High (>50,000) 3.115385 1.0000
## Atresia-Low (~1,000s) 2.740385 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# K+
# Linear Model
potassium_amb_fec_no_atresia_lm <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(potassium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.142711 -0.101391 -0.008131 0.034327 0.286925
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 2.6396171 0.0943482 27.977
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0004978 0.0005374 0.926
## Pr(>|t|)
## (Intercept) 1.38e-07 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.39
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1473 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.1251, Adjusted R-squared: -0.02069
## F-statistic: 0.8581 on 1 and 6 DF, p-value: 0.39
potassium_ambient_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm)
plot(potassium_ambient_sim)
## qu = 0.25, log(sigma) = -4.308757 : outer Newton did not converge fully.
## qu = 0.75, log(sigma) = -3.947908 : outer Newton did not converge fully.
plotResiduals(potassium_amb_fec_no_atresia_lm)
## qu = 0.25, log(sigma) = -4.308757 : outer Newton did not converge fully.
## qu = 0.75, log(sigma) = -3.947908 : outer Newton did not converge fully.
testDispersion(potassium_amb_fec_no_atresia_lm)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(chloride_ambient_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: chloride_ambient_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(potassium_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 2 3.612478 0.01534 0.12272
# Shapiro Test
shapiro.test(residuals(potassium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(potassium_amb_fec_no_atresia_lm)
## W = 0.87918, p-value = 0.1849
ggqqplot(residuals(lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root
potassium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(K...mmol.L. = sqrt(K...mmol.L.))
# Rerun lm model with square root transformed data
potassium_amb_fec_no_atresia_lm_sqrt <- lm(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
summary(potassium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = potassium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.043790 -0.030608 -0.001939 0.010984 0.085420
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.6239747 0.0283420 57.299
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0001535 0.0001614 0.951
## Pr(>|t|)
## (Intercept) 1.9e-09 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.378
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04424 on 6 degrees of freedom
## Multiple R-squared: 0.131, Adjusted R-squared: -0.01382
## F-statistic: 0.9046 on 1 and 6 DF, p-value: 0.3783
potassium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = potassium_amb_fec_no_atresia_lm_sqrt)
plot(potassium_amb_fec_sqrt_sim)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.
plotResiduals(potassium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.25. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
##
## DHARMa: qgam was unable to calculate quantile regression for quantile 0.5. Possibly to few (unique) data points / predictions. The quantile will be ommited in plots and significance calculations.
## qu = 0.75, log(sigma) = -3.903504 : outer Newton did not converge fully.
testDispersion(potassium_amb_fec_no_atresia_lm_sqrt)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
potassium_amb_fec_no_atresia_lm_res <- residuals(potassium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(potassium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: potassium_amb_fec_no_atresia_lm_res
## W = 0.88746, p-value = 0.2216
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(K...mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = potassium_ambient_sqrt)
Residuals Plot (if significant)
# K+
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
potassium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Weight Adjusted Fecundity", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(potassium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Calcium Graphs
# Ca+2 boxplot
calcium_ambient_fec_class_boxplot <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.), fill = "#189bff") +
geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Fecundity", y = "Calcium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
title = element_text(size = 12),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 10),
axis.text.y = element_text(size = 20),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 20))
print(calcium_ambient_fec_class_boxplot)
# Calcium Ambient No Atresia scatterplot
calcium_amb_fec_no_atresia_scatterplot <- Ambient_Samples %>%
ggplot() +
geom_point(aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.), shape = 21, fill = "#189bff") +
labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)") +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 30, hjust = 0.5),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_amb_fec_no_atresia_scatterplot)
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Differences: Between Fecundity Class
# Ca+2
# Ambient ANOVA Test
# Testing for differences in ambient parturition success metrics
# One-way ANOVA (parametric) or Kruskal-Wallis test (nonparametric)
calcium_aov_ambient_fec <- aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
summary(calcium_aov_ambient_fec)
## Df Sum Sq Mean Sq F value Pr(>F)
## Fecundity_Class 2 0.1328 0.06641 7.281 0.00482 **
## Residuals 18 0.1642 0.00912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# If significant, Post Hoc Test:
# - Tukey's HSD (parametric, normality & equal variance across groups)
# - Dunn's test (non-parametric, after Kruskal-Wallis test)
TukeyHSD(calcium_aov_ambient_fec)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## $Fecundity_Class
## diff lwr upr p adj
## Low (~1,000s)-High (>50,000) -0.0375000 -0.2098457 0.134845727 0.8450903
## Atresia-High (>50,000) -0.1807692 -0.3201293 -0.041409176 0.0103340
## Atresia-Low (~1,000s) -0.1432692 -0.2826293 -0.003909176 0.0433465
ca_ambient_fec_sim <- simulateResiduals(fittedModel = calcium_aov_ambient_fec)
## Warning in checkModel(fittedModel): DHARMa: fittedModel not in class of
## supported models. Absolutely no guarantee that this will work!
plot(ca_ambient_fec_sim)
plotResiduals(ca_ambient_fec_sim)
testDispersion(ca_ambient_fec_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91709, p-value = 0.872
## alternative hypothesis: two.sided
testOutliers(ca_ambient_fec_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: ca_ambient_fec_sim
## outliers at both margin(s) = 0, observations = 21, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.1610976
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(calcium_aov_ambient_fec)
## rstudent unadjusted p-value Bonferroni p
## 14 3.719954 0.0017024 0.035751
# Parametric tests
# Normality
shapiro.test(Ambient_Samples$Ca....mmol.L.)
##
## Shapiro-Wilk normality test
##
## data: Ambient_Samples$Ca....mmol.L.
## W = 0.9463, p-value = 0.2894
shapiro.test(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)))
##
## Shapiro-Wilk normality test
##
## data: residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))
## W = 0.91863, p-value = 0.08147
# QQplot:
calcium_ambient_fec_res_qqplot <- ggqqplot(residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))) +
labs(title = "Ambient Calcium Residual QQplot",
subtitle = "residuals(aov(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples))",
x = "Calcium Theoretical Quantiles (Predicted values)",
y = "Calcium Sample Quantiles") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
print(calcium_ambient_fec_res_qqplot)
# Parametric Variance test (Bartlett's test): require at least 2 obs for each group
bartlett.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Bartlett test of homogeneity of variances
##
## data: Ca....mmol.L. by Fecundity_Class
## Bartlett's K-squared = 6.1296, df = 2, p-value = 0.04666
# Nonparametric variance test (Levene's test):
leveneTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 1.2704 0.3047
## 18
# Nonparametric stat test: Kruskall Wallis
kruskal.test(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Kruskal-Wallis rank sum test
##
## data: Ca....mmol.L. by Fecundity_Class
## Kruskal-Wallis chi-squared = 9.153, df = 2, p-value = 0.01029
# Nonparametric Post Hoc: Dunns test
DunnTest(Ca....mmol.L. ~ Fecundity_Class, data = Ambient_Samples)
##
## Dunn's test of multiple comparisons using rank sums : holm
##
## mean.rank.diff pval
## Low (~1,000s)-High (>50,000) -2.500000 0.5686
## Atresia-High (>50,000) -9.528846 0.0216 *
## Atresia-Low (~1,000s) -7.028846 0.0948 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlations: Parameter to Numerical Fecundity (without Atresia IDs)
# Ca+2
# Linear Model
calcium_amb_fec_no_atresia_lm <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)
summary(calcium_amb_fec_no_atresia_lm)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = Ambient_Samples)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.02236 -0.01883 -0.01495 0.01027 0.06475
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.3912642 0.0207348 67.098
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0002731 0.0001181 2.313
## Pr(>|t|)
## (Intercept) 7.37e-10 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0601 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.03237 on 6 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.4713, Adjusted R-squared: 0.3831
## F-statistic: 5.348 on 1 and 6 DF, p-value: 0.06006
calcium_ambient_lm_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm)
plot(calcium_ambient_lm_sim)
## qu = 0.25, log(sigma) = -3.402339 : outer Newton did not converge fully.
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## qu = 0.75, log(sigma) = -2.863552 : outer Newton did not converge fully.
plotResiduals(calcium_ambient_lm_sim)
## qu = 0.25, log(sigma) = -3.402339 : outer Newton did not converge fully.
## Warning in newton(lsp = lsp, X = G$X, y = G$y, Eb = G$Eb, UrS = G$UrS, L = G$L,
## : Fitting terminated with step failure - check results carefully
## qu = 0.75, log(sigma) = -2.863552 : outer Newton did not converge fully.
testDispersion(calcium_ambient_lm_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
testOutliers(calcium_ambient_lm_sim)
##
## DHARMa outlier test based on exact binomial test with approximate
## expectations
##
## data: calcium_ambient_lm_sim
## outliers at both margin(s) = 0, observations = 8, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.007968127
## 95 percent confidence interval:
## 0.0000000 0.3694166
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 )
## 0
outlierTest(calcium_amb_fec_no_atresia_lm)
## No Studentized residuals with Bonferroni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferroni p
## 3 4.032387 0.0099976 0.079981
# Shapiro Test
shapiro.test(residuals(calcium_amb_fec_no_atresia_lm))
##
## Shapiro-Wilk normality test
##
## data: residuals(calcium_amb_fec_no_atresia_lm)
## W = 0.77177, p-value = 0.01417
ggqqplot(residuals(lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = Ambient_Samples)))
# Transform data: square root
calcium_ambient_sqrt <- Amb_Fec_No_Atresia_Samples %>%
mutate(Ca....mmol.L. = sqrt(Ca....mmol.L.))
# Rerun lm model with square root transformed data
calcium_amb_fec_no_atresia_lm_sqrt <- lm(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
summary(calcium_amb_fec_no_atresia_lm_sqrt)
##
## Call:
## lm(formula = Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight,
## data = calcium_ambient_sqrt)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.009302 -0.007894 -0.006232 0.004396 0.026830
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 1.180e+00 8.621e-03 136.816
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 1.143e-04 4.911e-05 2.328
## Pr(>|t|)
## (Intercept) 1.03e-11 ***
## Weight_Adjusted_Fecundity_Fecundity.Mean_Weight 0.0588 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.01346 on 6 degrees of freedom
## Multiple R-squared: 0.4745, Adjusted R-squared: 0.3869
## F-statistic: 5.418 on 1 and 6 DF, p-value: 0.05882
calcium_amb_fec_sqrt_sim <- simulateResiduals(fittedModel = calcium_amb_fec_no_atresia_lm_sqrt)
plot(calcium_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.
plotResiduals(calcium_amb_fec_sqrt_sim)
## qu = 0.25, log(sigma) = -3.52446 : outer Newton did not converge fully.
testDispersion(calcium_amb_fec_sqrt_sim)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.84221, p-value = 0.92
## alternative hypothesis: two.sided
# Residuals of model for shapiro test
calcium_amb_fec_no_atresia_lm_res <- residuals(calcium_amb_fec_no_atresia_lm_sqrt)
# Normality test
shapiro.test(calcium_amb_fec_no_atresia_lm_res)
##
## Shapiro-Wilk normality test
##
## data: calcium_amb_fec_no_atresia_lm_res
## W = 0.77335, p-value = 0.01475
# Parametric variance test: bartlett's test
# Require at least 2 obs for each group
# bartlett.test(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
# Nonparametric variance test (Levene's test)
# NA for quantitative explanatory variables
# leveneTest(Ca....mmol.L. ~ Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, data = calcium_ambient_sqrt)
Residuals Plot (if significant)
# Ca+2
# Ambient samples
# ggplot with geom_smooth(aes(x=,y=)), method = "lm")
calcium_ambient_lm_scatterplot <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Weight Adjusted Fecundity", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"),
plot.title = element_text(size = 25, hjust = 0.5),
axis.text.x = element_text(size = 15),
axis.text.y = element_text(size = 20),
axis.title.x = element_text(size = 20),
axis.title.y = element_text(size = 20))
print(calcium_ambient_lm_scatterplot)
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
Combined boxplots
# Fecundity class boxplots
ph_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = pH)) +
geom_point(aes(x = Fecundity_Class, y = pH), alpha = 0.5, colour = "black") +
labs(title = "pH", x = "Fecundity Class", y = "pH") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
hct_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Hct....)) +
geom_point(aes(x = Fecundity_Class, y = Hct....), alpha = 0.5, colour = "black") +
labs(title = "Hematocrit", x = "Fecundity Class", y = "Hematocrit (%)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
glu_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Glu..mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Glu..mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Glucose", x = "Fecundity Class", y = "Glucose (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
na_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Na...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Na...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Sodium", x = "Fecundity Class", y = "Sodium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
cl_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Cl...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Cl...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Chloride", x = "Fecundity Class", y = "Chloride (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
k_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = K...mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = K...mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Potassium", x = "Fecundity Class", y = "Potassium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
ca_fec_class <- ggplot(data = Ambient_Samples) +
geom_boxplot(aes(x = Fecundity_Class, y = Ca....mmol.L.)) +
geom_point(aes(x = Fecundity_Class, y = Ca....mmol.L.), alpha = 0.5, colour = "black") +
labs(title = "Calcium", x = "Fecundity Class", y = "Calcium (mmol/L)") +
facet_wrap(~ Ambient.Or.OAH) +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
fec_class_boxplots <- ph_fec_class + hct_fec_class + glu_fec_class + na_fec_class + cl_fec_class + k_fec_class + ca_fec_class + plot_layout(guides = "collect")
fec_class_boxplots
All Scatterplots
# lm plots
ph_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = pH)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "pH", x = "Weight Adjusted Fecudnity", y = "pH", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
hct_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Hct....)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Hematocrit", x = "Weight Adjusted Fecudnity", y = "Hematocrit (%)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
glucose_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Glu..mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Glucose", x = "Weight Adjusted Fecudnity", y = "Glucose (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
sodium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Na...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Sodium", x = "Weight Adjusted Fecudnity", y = "Sodium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
chloride_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Cl...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Chloride", x = "Weight Adjusted Fecudnity", y = "Chloride (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
potassium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = K...mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Potassium", x = "Weight Adjusted Fecudnity", y = "Potassium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
calcium_am_lm <- ggplot(data = Ambient_Samples, aes(x = Weight_Adjusted_Fecundity_Fecundity.Mean_Weight, y = Ca....mmol.L.)) +
geom_point() +
geom_smooth(aes(group = 1), method = "lm") +
labs(title = "Calcium", x = "Weight Adjusted Fecudnity", y = "Calcium (mmol/L)", color = "Category") +
scale_y_continuous() +
theme_classic() +
theme(panel.background = element_rect(fill = "white"),
plot.background = element_rect(fill = "white"),
legend.background = element_rect(fill = "white", color = "black"))
ambient_lm_plots <- ph_am_lm + hct_am_lm + glucose_am_lm + sodium_am_lm + chloride_am_lm + potassium_am_lm + calcium_am_lm + plot_layout(guides = "collect")
ambient_lm_plots
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 13 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Removed 13 rows containing missing values or values outside the scale range
## (`geom_point()`).